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UCTION 

The  parabolic  equation’  (PE)  method  requires  that  the 
acoustic  field  be  specified  at  some  range.  The  Gaussian  PE 
starter’  ’  and  more  accurate  PE  starters'"’  have  been  used 
widely  to  initialize  the  PE.  These  approaches  are  sufficient 
for  solving  two-dimensional  time-harmonic  propagation 
problems,  but  they  are  difficult  to  apply  to  three-dimension¬ 
al,  time-domain,  and  scattering  problems.  In  this  article,  we 
consider  a  PE  starter  that  is  accurate,  easy  to  construct  and 
understand  physically,  and  useful  for  both  frequency-do¬ 
main  and  time-domain  problems  in  both  two  and  three  spa¬ 
tial  dimensions  including  scattering  problems. 

Refraction  in  the  ocean  is  weak  and  thus  can  be  neglect¬ 
ed  near  a  sound  source  for  a  fairly  large  distance.  Reflections 
from  the  ocean  bottom  can  also  be  neglected  near  the  source 
because  the  reflection  coefficient  at  normal  incidence  is 
small.  Reflection  from  the  ocean  surface,  however,  cannot  be 
neglected  because  total  reflection  is  usually  assumed  at  the 
ocean  surface.  Thus  the  homogeneous  half-space  field  is  a 
valid  approximation  near  a  sound  source.  Far  from  the 
source,  both  refraction  and  reflection  from  the  ocean  bottom 
must  be  accounted  for.  Thus  the  half-space  field  breaks 
down  in  the  farfield,  and  a  boundary  layer’  exists  near  the 
source. 

We  use  the  method  of  stationary  phase”  to  show  that  the 
half-space  field  projects  properly  onto  the  normal  modes 
near  the  source.  Thus  the  half-space  field  can  be  used  to 
replace  a  source  condition  with  a  boundary  condition  on  a 
cylinder  enclosing  the  source  to  obtain  an  accurate  farfield 
solution.  We  compare  the  half-space  starter  and  the  Gaus¬ 
sian  PE  starter  in  both  the  frequency  domain  and  the  time 
domain.  For  deep-water  problems,  the  boundary  layer  be¬ 
havior  motivates  an  approach  for  speeding  up  PE  calcula¬ 
tions  as  well  as  an  efficient  ray-tracing  model. 

I.  THE  NORMAL-MODE  SOLUTION 

A  time-harmonic  steady  state  is  assumed,  and  the 
acoustic  pressure/?  is  factored  as /?(x,/)  =P(x)exp(  —/&>/), 
where  t  is  time,  x  is  the  Cartesian  position  v  ccior,  and  co  is  the 
circular  frequency.  The  complex  pressure  P  is  assumed  to 


satisfy  the  pressure  release  boundary  condition_.P  =  0  at  the 
ocean  surface,  the  outgoing  radiation  condition  at  infinity, 
and  the  reduced  wave  equation’ 

pV-[  ( 1//))VP] -f  A' -’P=  -  4nS(\  -  \i,),  (1) 

where  the  point  x„  is  the  source  location.  The  complex  wave- 
number  K  —  k(  I  +  icr/3)  is  used  to  model  sediment  loss, 
where  the  wavenumber  is  /c  =  a)/c,  a  —  (40m  log,,;  e)~\  P 
is  the  attenuation  in  decibels  per  wavelength  (dB/A),  and  c 
is  the  sound  speed.  The  variable  density  term  is  due  to  Berg- 
mann.” 

Since  the  ocean  is  a  waveguide,  energy  propagating  from 
a  source  exhibits  cylindrical  spreading  in  the  absence  of  azi-  ' 
muthal  variation  in  the  ocean.  Thus  it  is  often  beneficial  to 
solve  Eq.  ( 1 )  in  cylindrical  coordinates,  with  z  being  the 
depth  below  the  ocean  surface  and  r  being  the  horizontal 
distance  from  a  source  at  the  depth  z^y  Variations  in  both 
azimuth  and  range  are  assumed  to  be  sufficiently  weak  so 
that  the  region  near  the  source  can  be  treated  as  stratified, 
which  simplifies  Eq.  ( 1 )  to 

d-P  \  dp  dP  ^  d-P  ,  1  I 

dz~  p  dz  dz  (9m  r  dr 

=  -  —  3(r)S(z  -  z„).  (2) 

r 

The  normal-mode  solution’  of  Eq.  (2),  which  is  valid  in 
the  limit  1  and  1,  is 

,  imj  ^  ,  exp[((A:„  -  y,,)/-] 

P(r,z)~  I -  y  ((l„(z„)^„(z) - -  -  - 

V  '■  {K, 

(3) 

Kdz"  p  dz  dz  / 

where  (#„  (0)  =  (^„  ( oo  )  =0,  and  N  is  the  number  of  modes 
in  the  discrete  spectrum  of  L.  The  normal  modes  satisfy 
i<l>m  ^4>„ )  =  d,„„ ,  where  the  inner  product  associated  with  L 
is  defined  by 


Jn  p{z) 

The  constants  y„  are  defined  by 
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y„  =  k^,a{[3.(b],),  (6) 

where  A',,  =  A(z„).  The  vertical  propagation  angles  a„  are 
defined  by 

A-,,  =  A„cosa„.  (7) 

In  the  Appendix,  we  present  a  direct  asymptotic  derivation 
of  Eq.  (6),  which  has  also  been  derived  using  the  Rayleigh 
quotient.’ 

We  do  not  consider  frequencies  less  than  w,„  the  fre¬ 
quency  below  which  L  does  not  have  a  discrete  spectrum. 
For  each  frequency  oj  >  w„.  we  define  /I  (ru)  to  be  the  maxi¬ 
mum  of  a„  over  n.  We  define  e  to  be  the  maximum  of 
tan'  AUo)  over  (o.  The  discontinuity  in  sound  speed  at  the 
ocean  bottom  z  =  is  relatively  small.  Thus  energy  prop¬ 
agation  is  limited  to  small  vertical  angles,  and  we  assume 
that  1. 


II.  THE  HOMOGENEOUS  HALF-SPACE  FIELD 

Equation  (2)  is  difficult  to  solve  numerically  near  the 
source  because  P  and  the  Laplacian  are  singular  at  r  =  0. 
This  problem  can  be  handled  by  solving  Eq.  (2)  asymptoti¬ 
cally  near  the  source  for  r<r„.  Any  function  that  properly 
excites  the  normal  modes  can  be  used  as  a  boundary  condi¬ 
tion  at  r  =  r,,  to  replace  the  source  condition  at  r  =  0.  In  the 
high-frequency  limit,  it  is  necessary  to  accurately  account 
only  for  the  trapped  rays  that  propagate  within  the  angle /I  of 
the  horizontal;  and  the  homogeneous  half-space  field  pro¬ 
jects  properly  onto  the  normal  modes  at  r  =  r„  if  trapped 
rays  do  not  reflect  from  the  ocean  bottom  and  are  not  signifi¬ 
cantly  affected  by  refraction  in  the  ocean  for  r<r„.  This 
asymptotic  solution  breaks  down  for  large  r,,.  Thus  a  bound¬ 
ary  layer  exists  near  r  =  0,  and  the  half-space  solution  is  an 
inner  solution.  We  will  use  the  method  of  stationary  phase  to 
demonstrate  the  validity  of  the  half-space  solution  and  to 
show  when  it  is  valid. 

In  the  homogeneous  half-space  withy9  =  0,  the  solution 
of  Eq.  (2)  is 

P^(r,z)  =  e\p{ikd  _  )/d_  —  e\p{ikd  ^  )/d+  ,  (8) 

d"\  =  r  +  {z±z^y-.  (9) 


Since  a  ray  propagating  with  small  vertical  angle  a  remains 
'\nz<d  for  small  r,  one  would  expect  that 

(P,„d>J~iPA)-  (10) 


Since /9(z)  =  0  forz<ciand  <{>„  (z)  is  evanescent  forz>  tf,  it 
follows  from  Eq.  (6)  that  A„  'y,,  4,aP^\.  We  assume 
y„r4 1  <k„r,  and  the  factor  exp(  —  y„r)  may  be  dropped  in 
Eq.  (3).  We  obtain 


{P,d>„)- 


^  kr  Cos  a„ 


)exp  (/Arcosa,,).  (II) 


We  establish  Eq.  ( 10)  for  the  case  in  which  c  is  constant 
in  the  water  column.  Thus 


<^„(z)  =a,  sin  tAzsina,,)  (12) 

for  z<d,  where  a„  is  a  normalization  constant.  Because 
<l>„(z)  is  evanescent  for  z>d, 
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{P,,.6„)  (a„/2/)[/,  (r)  -  I  (r)  ~  J  ^(r)  +  J  (r)], 

(13) 


/.  (r) 

J  .  (r) 


exp{/A  [  \  r  -h  (z  —  z,,)’  +  z  sin  a„  ]  } 
+  (z-  z,,)- 

exp{/A  [  v7^  -f-  (z  -I-  z,,)^  ±  z  sin  a„  ] } 
\'r  +  (z 


dz , 

(14) 

dz . 

(15) 


Since  the  first  mode  is  essentially  half  of  a  cycle  of  a 
sinusoid  for  z<d,  we  deduce  from  Eq.  (12)  that 
Arf  sin  a,  =  0(  1 ).  This  implies  that  kd^l  because 
sin  1.  Thus/ .  and  7.  may  be  evaluated  by  the  method 
of  stationary  phase  with  the  length  scale  taken  to  be  d.  First, 
we  consider 

/_  =  r  B{z)e\p[ikil’(z)]dz,  (16) 

Jo 

il>=  \  r  -I-  (z  —  z,,)^  -  z  sin  a„,  (17) 


lA'  =  ^  -  - — sina„,  (18) 

fr  -f-  (z  -z,j)^ 

lA"  =  r[r -I- (z  -  z,,)-]  (19) 

From  Eq.  (18),  tA'  vanishes  at  z,  =  z„  4-  r  tan  a„ .  This  is  a 
stationary  point  of I_  if z,  <.d,  which  occurs  if 

r<r„  =  (d  -  Zt,)cota„.  (20) 

This  stationary  point  corresponds  to  the  ray  propagating 
downward  at  the  angle  a„.  At  the  stationary  point, 

B  =  r~'  cos 0!„  ,  (21 ) 

tA  =  r  cos  a„  —  z„  sin  a„  ,  (22) 


tA"  =  (r  cos'  a„)  '  .  (23) 

Thus  we  obtain 

/_  Im'/kr  cos  a„ 

Xexp(/Arcosa„  — /Az,,  sin  a„ ).  (24) 

A  similar  analysis  for  /+  gives  z,  =  z„  —  r  tan  a„ .  If 
z,  >0,  this  is  a  stationary  point  of  If  z,  <0, 1+  does  not 
have  a  stationary  point,  and  has  the  stationary  point  —  z, 
=  —z„  +  r  tan  a„ .  These  stationary  points  correspond  to 
the  ray  propagating  upward  at  the  angle  a„ ;/+(/+)  has  the 
stationary  point  if  the  ray  has  not  (has)  reflected  from  the 
ocean  surface.  If  z,  =  0,  /+  and  J+  share  this  stationary 
point.  Here,  J_  never  has  a  stationary  point.  In  each  case,  we 
find  that 


/+.(r)  —J+{r)  +J_{r) 

~yj2iri/kr  cosa„ 

Xexp(/Ar  cos  a„  -|- /Az,,  sin  a„ ),  (25) 

{P^,<i>„)~yj  liri/kr  cos  a„ 

Xa„  sin(Az„  sin  a„  )exp(/Arcos  a„).  (26) 

Thus  P/,  properly  excites  the  normal  modes. 

The  above  analysis  is  easily  modified  to  handle  the  case 
in  which  c  varies  in  the  water;  however,  several  cases  must  be 
considered.  If  <A„  (z)  is  not  evanescent  near  z  =  z,„  then  a„  is 
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real,  and  the  following  expansion  is  valid  forZ'-r,,; 

(2)  --  b„  s\n{k,^  sin  a„  +  d„)  .  (27) 

where  6.,  andr/„  are  constants  with  Ov,  d„  <  2 n-.  The  station¬ 
ary  phase  analysis  goes  through  similarly  in  this  case.  How¬ 
ever,  r  must  be  chosen  sufficiently  small  so  that  Eq.  (27)  is 
valid  at  the  stationary  points.  If  6„(z)  is  evanescent  near 
2  =  2,|,  it  suffices  to  show  that  the  projection  of  /*,,  onto  <!>„  is 
subdominant  to  the  projection  of  P,,  onto  other  modes;  d)„  is 
not  incorporated  into  the  exponentials  in  Eqs.  (14)  and  (15) 
in  this  case.  A  similar  analysis  results  in 


{P„.(t)„)  -  \  liri/krd),,  (z„)e\pUkr).  (28) 


Since  <A„  (2,1)  1,<E,,,(^„ )  is  subdominant.  It  is  possible  that 

Eq.  (28)  is  not  valid  if  the  main  contribution  of  one  of  the 
integrals  does  not  come  from  a  stationary  point.  However, 
the  contribution  is  automatically  subdominant  in  this  case. 

According  to  Eq.  (20),  r  must  be  small  enough  so  that  a 
ray  from  the  source  propagating  at  the  angle  a„  does  not 
intersect  the  sediment.  In  the  geometrical  optics  limit,  is 
the  range  at  which  P,,  breaks  down  for  rays  propagating  with 
vertical  angles  a  <  a„ .  For  lower  frequencies,  one  would  ex¬ 
pect  Pi,  to  break  down  for  smaller  r.  The  stationary  phase 
approximation  is  valid  if  exp(ikt//)  oscillates  rapidly  away 
from  the  stationary  point.  This  occurs  if 


=k 


2-2() 


—  sin  a„ 


,r  4-  (z-z,,)- 


(29) 


for  z  —  zj  =0(1).  Since  t/f  (z)  =  0(sin  a„ )  or  larger  for 
z  -  z, !  =  0(  1 ),  Eq.  (29)  holds  for  the  case  k  sin  a„  >  1. 
For  k  sin  a„  =0(1),  Eq.  (29)  is  equivalent  to 


A- iz-ZoiAr -4  (2-z„)->l.  (30) 

Let  r  =  Kr„  =  K(d  —  Z||)cot  a„,  where 

A-(c(-z„)<^l.  (31) 

From  Eq.  (31),  (z,  —  z„|  = /rid  —  z,,! 1.  We  deduce  that 
z  —  z,  I  =  0(  1 )  is  equivalent  to  |z  —  zJ  =0(1),  and  the 
stationary  phase  approximation  is  valid  ifEq.  (30)  holds  for 
'z  —  z,i)  =0(1).  For  r  =  0(  1 ),  Eq.  (30)  holds  since  A>  1. 
Due  to  the  fact  that  A  sin  a„  =0(1),  Eq.  (30)  follows  from 
Eq.  ( 31 )  for  r>  1.  Thus  the  stationary  phase  evaluation  of 
the  integrals  in  Eqs.  (14)  and  (15)  is  valid  for  r4r„  if 
d  —  z,i  =  Oi\}  and  for  r<r„  if  d  —  z,,  <  1  or  A  sin  a„  >  1 . 


III.  APPLICATIONS 

With  the  inner  solution,  it  is  possible  to  avoid  the  singu¬ 
lar  nature  of  Eq.  (2)  near  r  =  0.  The  source  condition  at 
r  =  0  is  replaced  by  the  boundary  condition  P(r„,z) 
=  (r,„z).  There  is  an  additional  problem  to  be  dealt  with, 

however,  before  the  outer  problem  can  be  solved  numerical¬ 
ly  in  the  outer  region  r>  r„.  The  infinite  spatial  domain  must 
be  truncated.  The  outer  solution  of  Eq.  (2)  replaces  the 
boundary  value  problem  with  an  initial  value  problem  in 
range  taking  care  of  half  of  this  problem.  The  problem  is 
truncated  in  depth  with  the  boundary  condition 
P{r,z^f )  =  0,  where  Zy,  is  large.  Since  the  sediment  is  lossy, 
only  a  negligible  amount  of  energy  can  travel  to  the  lower 
boundary  and  return  to  the  water  column.  Thus  only  a  small 
error  is  introduced  by  truncating  the  domain.  The  inner  so¬ 


lution  P/,  does  not  satisfy  the  pressure  release  boundary  con¬ 
dition  at  2  =  Zv,.  resulting  in  Gibbs’  oscillations.  However, 
the  propagating  modes  in  the  truncated  dumain  aic  evanes¬ 
cent  for  z:>d.  Thus  the  Gibbs'  oscillations  project  onto  non¬ 
propagating  modes  and  decay  rapidly  with  r  causing  little 
error  in  the  farfield. 

The  major  benefit  of  the  half-space  field  is  that  it  can  be 
applied  to  several  problems.  If  azimuthal  variations  are  as¬ 
sumed  to  be  a  weak  perturbation,  the  half-space  field  can  be 
applied  to  initialize  the  three-dimensional  parabolic  equa¬ 
tion  at  r  =  /■„.  For  a  pulsed  source  with  the  source  function 
/(/)  of  finite  bandwidth,  the  homogeneous  half-space  field 


p„(r,z.t) 

d  \  c„  1 


(32) 


is  '’?!id  inner  sohi*i'sn  of  the  time-domain  problem  by  the 
superposition  principle  since  it  is  valid  for  each  frequency 
composing  the  source.  The  source  condition  at  /•=  0  is  re¬ 
placed  by  the  boundary  condition  /)(r,„z,/)  =  p,,  (r,|,z,t). 
This  starter  has  been  applied  to  initialize  the  time-domain 
parabolic  equation’*  (TDPE)  and  the  progressive  wave 
equation'*’  producing  accurate  results. 

The  inner  solution  can  also  be  applied  to  simplify  the 
waveguide  scattering  problem  for  a  bounded  object  in 
Q)<z<d.  Previous  work  on  this  problem  has  not  exploited 
asymptotic  limits."  ’’  The  total  field  P,  is  defined  to  be  the 
sum  of  the  incident  and  scattered  fields  P,  and  P,,  z„  is  as¬ 
sumed  to  be  a  parameter,  and  the  solution  of  Eq.  ( 2 )  is  writ¬ 
ten  as  P(r,z;Z|,).  The  scattered  field  has  the  representation 


C  r  ( dG  dP,  \ 

P-W  =  — (x,Xo)P,(x„)  -  — (x„)G(x,x„)  W, 

J  J  \  ufi  dfi  / 


(33) 


where  G  is  the  waveguide  Green’s  function  defined  by 

G(x,x„)  =  ( l/47r)P  [/U*-  X,,)-  -t-  (p  -p„)%z;Z||]  .  (34) 

The  integral  is  over  the  surface  of  the  scatterer.  Since  G  is  the 
field  due  to  a  point  source,  and  dG /dn  can  be  approximated 
by  two  point  sources.  P  may  be  regarded  as  the  field  due  to  a 
collection  of  point  sources  on  and  near  the  surface  of  the 
scatterer.  Thus  the  scatterer  behaves  like  a  collection  of 
point  sources,  which  suggests  that  a  valid  inner  solution  for 
the  scattering  problem  consists  of  the  solution  of  the  scatter¬ 
ing  problem  in  the  homogeneous  half-space.  This  is  a  signifi¬ 
cant  simplification  of  the  waveguide  scattering  problem  be¬ 
cause  useful  methods  for  solving  the  half-space  problem 
have  been  developed.  An  asymptotic  analysis  of  the  wave¬ 
guide  scattering  problem  is  given  in  Ref  13. 


IV.  COMPARISON  WITH  THE  GAUSSIAN  PE  STARTER 


There  is  an  alternate  method  for  dealing  with  the  singu¬ 
lar  behavior  near  the  source.  From  Eq.  (3),  it  follows  that 
P~r  ''-Q  in  the  farfield,  where 


d-Q  \  dp  dQ 
dz~  p  dz  dz 

Q{Q,z)  =  ^2iri  2 
»  1 


(35) 

(36) 
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To  av  oid  computing  normal  modes,  the  Gaussian  PE  starter 


-  A  ,^(z  +  z„)-j 


(37) 


which  approximates  £)(z,z„:Aii)  ~  (2~//A„) '  ’[6(z  —  z„) 
-  t>(z  +  Z|,)  ],  is  used  in  place  of  Eq.  (36)  as  an  initial  condi¬ 
tion  at  r  =  0.  This  is  a  valid  approximation  since  A„  ~  A,,  and 


{Q.Sj  =  ^2-i/k„<t>„iz„),  (38) 

{D.d,„)  =^2-i/k,,6„iz„).  (39) 


Since  k„  diverges  from  A„  as  a„  increases,  we  see  from 
Eqs.  (26)  and  (39)  that  P,,  excites  the  higher  modes  better 
than  D.  which  should  excite  modes  better  than  the  approxi¬ 
mation  r.  We  illustrate  this  with  an  example  for  which  data 
appear  in  Table  I.  The  subscript  w  stands  for  water  value. 
The  subscript  b  stands  for  bottom  value.  Projections  of  P/,, 
D.  and  T  onto  the  modes  appear  in  Table  II.  We  observe  that 
the  error  in  P,,  is  small  and  uniform  in  a„,  while  the  error  in 
r  increases  with  a„  and  is  larger  than  the  error  in  D  for  all 
a„. 

A  pulsed  source  function  is  decomposed  as 

fU)  =  f6i)exp(  —  icot)d(o. 

By  superposition,  p~r  '  in  the  farficld,  where 
d  ~q  1  dp  dq  ^  d'q  _  J_  d'q 

dz-  p  dz  dz  dr  r  dt' 

qiO.z.t)  =  j  F(cj)Q^.z\ —'j  e\p{  —  i(ot)dci),  (42) 

and  C|,  =  c(Z|,).  We  consider  the  source  function  /(/) 
=  exp  [  —  ( tv )  ■  I ,  for  which 

f(w)  =  (2i\-)  '  exp[  -  (w/2v)-l.  (43) 


(40) 

(41) 


Substituting  D  for  Q  in  Eq.  (42),  we  obtain 
q(0,zj)  =  —  _  [5(z  -  z,,)  -  S(z  +  Z|,)  ] 


exp(  —  icot)d(o. 

(44) 


The  delta  function  in  Eq.  (44)  can  be  approximated  by  a 
Gaussian  of  width  w: 


6(z)^(u\Tr)  '  exp[  -  (z/w)-\.  (45) 


Substituting  Eq.  (45)  in  Eq.  (44)  and  manipulating  the  inte¬ 
gral,  we  obtain  the  Gaussian  TDPE  starter: 


^(0,z./)  =  — 
wv 


jexpf  - 

-)1 

-expf 

1  L  \ 

1 

L 

f-i-exp 

2"  vw 

cos| 

(46) 


We  consider  the  example  described  in  Table  III.  The 
numerical  solution  of  the  TDPE  requires  that  q(r,z,l)—0 
rapidly  as/->  —  cc  (see  Ref.  9).  We  see  from  Fig.  1  that  this 
condition  holds.  A  sequence  of  contour  plots  generated  us¬ 
ing  the  TDPE  with  the  Gaussian  starter  appears  in  Fig.  2. 
From  the  plots  generated  using  the  half-space  starter,  which 
appear  in  Fig.  3,  we  deduce  that  the  Gaussian  starter  evolves 
properly.  The  accuracy  of  the  Gaussian  starter  is  illustrated 
quantitatively  in  Fig.  4.  The  error  in  the  surface-reflected 
arrival  near  t  —  60  ms,  which  propagates  at  roughly  1 5°,  is  on 
the  order  of  10%.  We  see  from  Table  II  that  this  error  is  of 
the  appropriate  magnitude. 


TABLE  II.  Projections  of  P,,.  T.  and  D  onto  modes. 


u 

0^ 

ID 

c/} 

(/> 

u 

y„  ( m  ' ) 

<P„.i.,) 

<P.<t'..) 

(DA.) 

(PA„> 

1.04676 

0.00000 

1.65951 

0.99901 

0.99895 

0.99919 

Qi 

a. 

1.04544 

0,00002 

3.31619 

1.01666 

0.99581 

0.99916 

Q  0- 

1.0432.V 

0.00004 

4.96723 

1.00280 

0,99060 

0.99812 

u 

L04014 

0.00(X)7 

6.60980 

0.98127 

0.98333 

0.99667 

1.0.1616 

0.00012 

8.24106 

0.99439 

0.97407 

0.99482 

1  0,1 1 26 

0.(XX)17 

9.85815 

1.02109 

0.96285 

0.99259 

tK 

1.02546 

0.00024 

11.45808 

1.00875 

0.94974 

0.98999 

Z 

1.01X74 

0.(XX).1,1 

13.03760 

0.97815 

0.9.1483 

0.98703 

1  01109 

0.00046 

14.59277 

0.98854 

0.91823 

0.98374 

300  C 

1.00252 

0.00068 

16.11757 

1.01759 

0.90008 

0.98015 

RELATIVE  TlME(ms) 

0.99314 

0.00126 

17.59651 

1.01015 

0.88069 

0.97632 

FIG.  1.  The  Gaussian  TDPE  starter  at  z  —  7.S  m. 
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I)F.PTH(in)  S  nKF'TH(iu)  3  2;  l)KlTll(ii>)  a;  DEPTH(m) 


j.  2.  TDPE  solution  generated  with  the  Gaussian  starter.  The  solid  contours  corresptrnd  \op>0.  The  dashed  contours  correspond  lop  <0.  (a)  r  =  0,  (b) 
50  m,  (c)  r  =  100  m,  (d)  r  =  500  m. 


FIG,  3.  TDPE  solution  generated  with  the  half-space  starter,  (a)  r=  50  m, 
(b)  r  =  KX)  m,  (c)  r  =  500  m.  The  agreement  with  the  solution  generated 
with  the  Gaussian  starter  is  qualitatively  good  at  r  =  500  m. 
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FIG-  4.  TDPE  solutions  at  r  —  5(X)  m  and  2  =  50  m  generated  with  the 
Gaussian  starter  (solid)  and  the  halFspace  starter  (dashed).  The 
Gaussiun-generaied  sidution  breaks  down  for  late  arrivals  that  propagate  at 
large  angles. 


V.  EFFICIENT  PE  AND  RAY-TRACING  CALCULATIONS 


The  PE  method  is  valid  for  high  frequencies.  However, 
the  computation  time  required  for  the  numerical  solution  of 
the  PE  is  proportional  to  the  frequency  squared  in  two  di¬ 
mensions.  To  avoid  time-consuming  calculations,  one  is 
tempted  to  apply  ray  theory.  However,  the  ray  theory  solu¬ 
tion  can  be  unpleasant  to  implement  due  to  the  fact  that  it  is 
implicit.  In  this  section,  we  discuss  approaches  for  speeding 
up  PE  calculations  and  for  obtaining  qualitative  results  effi¬ 
ciently  with  ray  tracing. 

One  of  the  standard  methods  for  matching  the  inner  and 
outer  solutions  is  to  obtain  the  behavior  of  the  inner  solution 
in  the  outer  region.'"'  The  half-space  field  P/,  exhibits  beams 
in  the  outer  region  due  to  the  Lloyd’s  mirror  effect.  Thus  one 
would  expect  the  outer  solution  to  consist  of  beams  that  are 
curved  due  to  'efraction.  We  illustrate  the  beams  in  the  outer 
region  for  a  200-Hz  source  placed  at  2,,  =  50  m  in  a  deep¬ 
water  ocean  with  Munk’s  exponential  sound-speed  profile'' 


1  +p 


H 
z  —  z, 


+  expi 


(47) 


where  /r^l,  z^.,,  is  the  channel  depth,  is  the  channel 
speed,  and  //is  the  channel  thickness.  A  plot  ofc„,  appears  in 
Fig.  5.  A  plot  of  PE-generated  transmission  loss  appears  in 


1575 


y. 


X  1500  — '  ' 

y 

1475  I - 1  '  I  -  -L 

0  1700  :)400  5100 

DEPTH(m) 

FIG.  5  Munk’s  deep-water  sound-speed  profile. 


0  ZO  40  60  80  100 


RANGE(km) 

FIG.  6.  PE-generated  transmission  loss.  The  curved  beams  of  the  outer  so¬ 
lution  can  be  traced  back  to  the  Lloyd's  mirror  beams  near  the  origin. 


Fig.  6  for  /U  =  0.0071,  c^.,,  =  1500  m/s,  z^^^  =  1000  m,  and 
H  =  1200  m.  The  curved  beams,  which  can  be  traced  to  the 
Lloyd’s  mirror  beams  near  the  source,  are  the  dominant  fea¬ 
ture.  The  beams  evolve  to  a  constant  width  in  the  two-di¬ 
mensional  vertical  cross  section,  which  is  necessary  since 
cylindrical  spreading  applies  in  the  sound  channel.  The  PE 
calculation  took  2.5  h  on  a  Digital  VAX-8650  computer. 

The  dipole  behavior  exhibited  in  Fig.  6  suggests  an  ap¬ 
proach  for  speeding  up  PE  calculations.  Since  the  beam  pat¬ 
tern  is  determined  by  transmission  loss  for  the  frequen¬ 
cy  (o  can  be  approximated  efficiently  by  replacing  co  with  fcj, 
z„  with  z,/f,  and  /?  with  /?/^,  where  f  <  1.  We  refer  to  this 
approach  as  the  f  method.  For  it  to  be  valid,  ^co  must  be  large 
enough  for  ray  theory  to  be  valid.  We  illustrate  the  ^method 
by  taking  f  =  1/2  for  the  2(X)-Hz  problem  with  Munk’s  pro¬ 
file.  The  transmission  loss  plot  obtained  by  the  ^  method 
appears  in  Fig.  7  and  agrees  well  with  Fig.  6.  Since  the  f 
method  gives  accurate  results  and  decreases  the  run  time  by 
the  factor  ^ "  in  two  dimensions  and  by  the  factor  f '  in  three 
dimensions,  it  should  be  useful  for  underwater  acoustics  cal¬ 
culations. 

Ray-tracing  calculations  are  complicated  because  it  is 
necessary  to  determine  all  of  the  rays  that  pass  through  a 
given  point  and  sum  their  contributions  to  the  field  at  that 
point.  A  simple  way  to  avoid  this  difficulty  is  to  take  into 
account  the  behavior  in  the  inner  region.  A  good  qualitative 
representation  of  the  farfield  is  obtained  by  tracing  rays  from 
the  origin  in  the  directions  of  the  Lloyd’s  mirror  beams,  as  in 
Fig.  8.  This  calculation  required  just  3  s  on  a  VAX-8650,  yet 
the  result  gives  a  useful  description  of  the  farfield. 


£ 


s 


FIG.  7.  The  C  method  with  C  =  1/2.  There  is  good  agreement  with  Fig.  6. 
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FIG.  8.  Rays  launched  in  the  directions  of  the  Lloyd’s  mirror  beams  to 
obtain  the  qualitative  behavior  in  Fig.  6. 

These  concepts  of  ray-tracing  and  matched  asymptotics 
combine  to  motivate  an  efficient  method  for  source  localiza¬ 
tion.  Suppose  that  a  high-frequency  (  >  100  Hz),  time-har¬ 
monic  signal  from  a  point  source  of  unknown  location  is 
received  at  a  vertical  array.  Signal  processing  techniques  can 
be  used  to  determine  the  depths  and  angles  of  beam  arrivals 
at  the  array.  To  determine  the  source  location,  rays  corre¬ 
sponding  to  the  beams  are  traced  away  from  the  array  in 
both  directions.  The  range  of  the  source  is  the  range  at  which 
all  of  the  rays  meet  at  one  point  on  the  ocean  surface.  The 
depth  of  the  source  is  determined  by  the  Lloyd’s  mirror 
beam  pattern  described  by  the  rays  at  this  range. 

To  further  illustrate  the  ray-tracing  method,  we  consid¬ 
er  an  ocean  with  the  Munk  profile  used  above  and  having  the 
bathymetry 


=d„-'^  h„  exp(  -  a;,)  , 

n 

(48) 

cr„  =  [(A -A,,)- 4-  (y-yJ-J/ic;  . 

(49) 

Data  for  this  problem  appear  in  Table  IV.  We  apply  the  ray¬ 
tracing  method  to  determine  the  locations  at  which  beams 
are  incident  upon  the  ocean  surface  at  the  first  two  conver¬ 
gence  zones  and  plot  them  in  Fig.  9.  Although  beams  lose 
energy  as  they  reflect  from  the  ocean  bottom,  they  retain 
some  energy,  as  is  evident  in  the  shadow  zone  forz  <  2500  m 
and  r  near  25  km  in  Fig.  6.  Thus  we  do  not  terminate  the  rays 
as  they  reflect  from  the  seamounts.  The  seamounts  result  in 
various  partial  and  total  shadow  zones.  Some  of  the  rays 
encounter  the  seamount  nearest  the  origin  two  times  produc¬ 
ing  additional  cusps  in  the  curves. 

VI.  CONCLUSIONS 

The  homogeneous  half-space  field  is  useful  for  time-har¬ 
monic  and  pulsed  point  and  multipole  sources  in  both  two 
and  three  spatial  dimensions.  It  is  probably  the  simplest  pos¬ 
sible  inner  solution  because  it  accounts  for  rays  from  both 
the  source  point  and  the  image  point  at  which  a  source  of 


TABLE  IV.  Data  for  the  three-dimensional  ray-trace  problem. 


X|  =  20  km 

X.  =  —  70  km 

j:,  =  10  km 

y,=0 

i»  =  ?0  km 

y.=  —  60  km 

A,  =  2500  m 

A,  =  3000  m 

A,  =  2500  m 

w,  =  10  km 

w.  =  15  km 

u),  =  20  km 

d„  =  5000  m 

(O  —  40077S  ‘ 

z„  =  50  m 

FIG.  9.  View  from  above  an  ocean  with  seamounts.  The  curves  mark  the  | 
locations  of  beam  arrivals  the  ocean  surface.  Various  shadow  zones  ap-  ,  ; 
pear  behind  the  seamounts.  '  j 

i 

opposite  phase  is  needed  to  account  for  the  pressure  release 
surface.  The  half-space  field  is  easy  to  interpret  physically. 
Refraction  is  weak  in  the  ocean  and  can  be  neglected  over  ■ 
short  ranges.  Reflections  from  the  ocean  bottom  occur  at 
large  angles  near  the  source  and  thus  do  not  affect  the  far- 
field.  These  concepts  lead  to  an  important  simplification  of 
the  scattering  problem.  The  concepts  of  matched  asymptot¬ 
ics  motivate  the  f  method  and  the  ray-tracing  model. 
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APPENDIX:  ASYMPTOTIC  DERIVATION  OF 
EIGENVALUES 

We  assume  the  following  ansatz,  which  is  similar  to  the 


WKBJ  ansatz,  for  the  eigenvalue  problem: 

P{r,z)-~m  2  d>Az„)<t>Az)H'„"{K„r),  (Al) 

n  rr  1 

<t>„(z)=fi}„(7)-f-i^„(z)-|-...,  (A2) 

—  ^11  +  ‘Yn  +  •  (A3) 

Substituting  the  «th  term  of  Eq.  (Al)intoEq.  (2),  we  ob¬ 
tain 

Lil>„  -  k  ltl>„  -  2/(  -  aPk  -  -|-  )d„.  (A4) 

Equation  ( A4)  has  a  solution  only  if  the  following  solvabil¬ 
ity  condition  is  satisfied: 

A  -  al3k- +  y„kn)d>, ,,<!>„)  =0,  (A5) 

Y„  =  ( a/k„ )  {13k  -,61 )  .  ( A6 ) 

Equation  (6)  follows  since  k~k,^  and  k„  -  k^y 
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